This file is used to analyse the IFE granular and spinous dataset.
library(dplyr)
library(patchwork)
library(ggplot2)
library(ComplexHeatmap)
.libPaths()
## [1] "/usr/local/lib/R/library"
In this section, we set the global settings of the analysis. We will store data there :
save_name = "ifegs"
out_dir = "."
We load the dataset :
sobj = readRDS(paste0(out_dir, "/", save_name, "_sobj.rds"))
sobj
## An object of class Seurat
## 15051 features across 994 samples within 1 assay
## Active assay: RNA (15051 features, 2000 variable features)
## 6 dimensional reductions calculated: RNA_pca, RNA_pca_22_tsne, RNA_pca_22_umap, harmony, harmony_22_umap, harmony_22_tsne
We load the sample information :
sample_info = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_sample_info.rds"))
project_names_oi = sample_info$project_name
graphics::pie(rep(1, nrow(sample_info)),
col = sample_info$color,
labels = sample_info$project_name)
Here are custom colors for each cell type :
color_markers = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_color_markers.rds"))
data.frame(cell_type = names(color_markers),
color = unlist(color_markers)) %>%
ggplot2::ggplot(., aes(x = cell_type, y = 0, fill = cell_type)) +
ggplot2::geom_point(pch = 21, size = 5) +
ggplot2::scale_fill_manual(values = unlist(color_markers), breaks = names(color_markers)) +
ggplot2::theme_classic() +
ggplot2::theme(legend.position = "none",
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 30, hjust = 1))
This is the projection of interest :
name2D = "harmony_22_tsne"
We design a custom function to make the GSEA plot and a word cloud graph :
make_gsea_plot = function(gsea_results, gs_oi, fold_change, metric = "FC") {
fold_change$metric = fold_change[, metric]
plot_list = lapply(gs_oi, FUN = function(gene_set) {
# Gene set content
gs_content = gene_sets %>%
dplyr::filter(gs_name == gene_set) %>%
dplyr::pull(ensembl_gene) %>%
unique()
# Gene set size
nb_genes = length(gs_content)
# Enrichment metrics
NES = gsea_results@result[gene_set, "NES"]
p.adjust = gsea_results@result[gene_set, "p.adjust"] %>%
round(., 4)
qvalues = gsea_results@result[gene_set, "qvalues"]
if (p.adjust > 0.05) {
p.adjust = paste0("<span style='color:red;'>", p.adjust, "</span>")
}
my_subtitle = paste0("\nNES : ", round(NES, 2),
" | padj : ", p.adjust,
" | qval : ", round(qvalues, 4),
" | set size : ", nb_genes, " genes")
# Size limits
lower_FC = min(fold_change[gs_content, ]$metric, na.rm = TRUE)
upper_FC = max(fold_change[gs_content, ]$metric, na.rm = TRUE)
# Plot
p = enrichplot::gseaplot2(x = gsea_results, geneSetID = gene_set) +
ggplot2::labs(title = gene_set,
subtitle = my_subtitle) +
ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold",
margin = ggplot2::margin(3, 3, 5, 3)),
plot.subtitle = ggtext::element_markdown(hjust = 0.5,
size = 10))
wc = ggplot2::ggplot(fold_change[gs_content, ],
aes(label = gene_name, size = abs(metric), color = metric)) +
ggwordcloud::geom_text_wordcloud_area(show.legend = TRUE) +
ggplot2::scale_color_gradient2(
name = metric,
low = aquarius::color_cnv[1],
mid = "gray70", midpoint = 0,
high = aquarius::color_cnv[3]) +
ggplot2::scale_size_area(max_size = 7) +
ggplot2::theme_minimal() +
ggplot2::guides(size = "none")
return(list(p, wc))
}) %>% unlist(., recursive = FALSE)
return(plot_list)
}
We visualize gene expression for some markers :
features = c("percent.mt", "nFeature_RNA", "TK1",
"SPINK5", "KRT1", "KRTDAP")
plot_list = lapply(features, FUN = function(one_gene) {
Seurat::FeaturePlot(sobj, features = one_gene,
reduction = name2D) +
ggplot2::theme(aspect.ratio = 1) +
ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
Seurat::NoAxes()
})
patchwork::wrap_plots(plot_list, ncol = 3)
We visualize clusters :
cluster_plot = Seurat::DimPlot(sobj, reduction = name2D, label = TRUE) +
Seurat::NoAxes() +
ggplot2::theme(aspect.ratio = 1)
cluster_plot
KRT1 and KRTDAP being specifically expressed in IFE spinous, we make a cluster annotation:
sobj$cluster_type = case_when(
sobj$seurat_clusters %in% c(5,6) ~ "IFE proliferative",
sobj$seurat_clusters %in% c(1,2,3) ~ "IFE spinous",
sobj$seurat_clusters %in% c(0,4) ~ "IFE granular",
TRUE ~ "others") %>% # others: 8,9,7,0
factor(., levels = c("IFE spinous", "IFE granular",
"IFE proliferative", "others"))
cluster_color = setNames(nm = c("IFE spinous", "IFE granular", "IFE proliferative", "others"),
c("steelblue4", "lightskyblue", "gray50", "gray10"))
We visualize the cluster annotation split by sample :
plot_list = aquarius::plot_split_dimred(sobj,
reduction = name2D,
split_by = "project_name",
group_by = "cluster_type",
split_color = setNames(sample_info$color,
nm = sample_info$project_name),
group_color = cluster_color,
bg_pt_size = 0.5, main_pt_size = 0.5)
plot_list[[length(plot_list) + 1]] = Seurat::DimPlot(sobj, reduction = name2D,
group.by = "cluster_type") +
ggplot2::scale_color_manual(values = cluster_color,
breaks = names(cluster_color)) +
Seurat::NoAxes() +
ggplot2::theme(aspect.ratio = 1)
patchwork::wrap_plots(plot_list, ncol = 4) &
Seurat::NoLegend()
We save the results in a list :
list_results = list()
We make over-representation analysis for each group of genes. We load gene sets from MSigDB :
gene_sets = aquarius::get_gene_sets(species = "Homo sapiens")
gene_sets = gene_sets$gene_sets
head(gene_sets)
## # A tibble: 6 x 16
## gs_cat gs_subcat gs_name gene_symbol entrez_gene ensembl_gene human_gene_symb~
## <chr> <chr> <chr> <chr> <int> <chr> <chr>
## 1 C5 GO:BP GOBP_1~ AASDHPPT 60496 ENSG0000014~ AASDHPPT
## 2 C5 GO:BP GOBP_1~ ALDH1L1 10840 ENSG0000014~ ALDH1L1
## 3 C5 GO:BP GOBP_1~ ALDH1L2 160428 ENSG0000013~ ALDH1L2
## 4 C5 GO:BP GOBP_1~ MTHFD1 4522 ENSG0000010~ MTHFD1
## 5 C5 GO:BP GOBP_1~ MTHFD1L 25902 ENSG0000012~ MTHFD1L
## 6 C5 GO:BP GOBP_1~ MTHFD2L 441024 ENSG0000016~ MTHFD2L
## # ... with 9 more variables: human_entrez_gene <int>, human_ensembl_gene <chr>,
## # gs_id <chr>, gs_pmid <chr>, gs_geoid <chr>, gs_exact_source <chr>,
## # gs_url <chr>, gs_description <chr>, category <chr>
How many gene sets ?
gene_sets[, c("gs_subcat", "gs_name")] %>%
dplyr::distinct() %>%
dplyr::pull(gs_subcat) %>%
table() %>%
as.data.frame.table() %>%
`colnames<-`(c("Category", "Nb gene sets"))
## Category Nb gene sets
## 1 50
## 2 CP:KEGG 186
## 3 CP:PID 196
## 4 CP:REACTOME 1615
## 5 CP:WIKIPATHWAYS 664
## 6 GO:BP 7658
## 7 GO:CC 1006
## 8 GO:MF 1738
We get gene name and gene ID correspondence :
gene_corresp = sobj@assays[["RNA"]]@meta.features[, c("gene_name", "Ensembl_ID")] %>%
`colnames<-`(c("NAME", "ID")) %>%
dplyr::mutate(ID = as.character(ID))
rownames(gene_corresp) = gene_corresp$ID
head(gene_corresp)
## NAME ID
## ENSG00000238009 AL627309.1 ENSG00000238009
## ENSG00000237491 AL669831.5 ENSG00000237491
## ENSG00000225880 LINC00115 ENSG00000225880
## ENSG00000230368 FAM41C ENSG00000230368
## ENSG00000230699 AL645608.3 ENSG00000230699
## ENSG00000187634 SAMD11 ENSG00000187634
We perform a differential expression analysis between granular and spinous clusters.
group_name = "granular_vs_spinous"
Seurat::Idents(sobj) = sobj$cluster_type
table(Seurat::Idents(sobj))
##
## IFE spinous IFE granular IFE proliferative others
## 350 260 153 231
We visualize these populations :
aquarius::plot_red_and_blue(sobj, reduction = name2D,
group1 = "IFE granular", group2 = "IFE spinous")
We identify DE genes between IFE granular and IFE spinous :
mark = Seurat::FindMarkers(sobj, ident.1 = "IFE granular", ident.2 = "IFE spinous")
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 132 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## PTN 2.884716e-60 1.2330983 0.927 0.340 4.341786e-56
## LGALS1 7.228537e-19 1.0315635 0.592 0.274 1.087967e-14
## APOE 5.473262e-42 0.7711590 0.992 0.886 8.237807e-38
## CHI3L1 7.302950e-33 0.7452335 0.631 0.160 1.099167e-28
## C1QTNF12 4.219449e-42 0.7192155 0.781 0.254 6.350692e-38
## CYP27A1 1.879473e-51 0.7034252 0.758 0.157 2.828795e-47
## DAPK2 3.697558e-49 0.7027778 0.808 0.246 5.565195e-45
## CLDN1 1.697175e-30 0.6445342 0.881 0.551 2.554419e-26
## EFNB2 9.527752e-43 0.6315371 0.765 0.206 1.434022e-38
## CYP1B1 3.413979e-35 0.6183131 0.565 0.091 5.138379e-31
## ADGRL3 1.717088e-43 0.5942463 0.788 0.226 2.584389e-39
## S100A2 1.920555e-18 0.5641180 0.992 0.931 2.890627e-14
## NFIB 2.544404e-23 0.5558578 0.923 0.657 3.829582e-19
## MALAT1 1.119576e-15 0.5345120 0.985 0.869 1.685074e-11
## CD109 1.279586e-32 0.5342700 0.835 0.460 1.925904e-28
## ZFP36L2 3.209914e-23 0.5297392 0.958 0.734 4.831241e-19
## ACTN1 2.516650e-25 0.4620441 0.881 0.589 3.787809e-21
## S100A6 2.860886e-18 0.4473052 0.977 0.869 4.305920e-14
## XIST 1.069519e-12 0.4310408 0.704 0.423 1.609734e-08
## APOC1 7.124088e-24 0.4264331 0.712 0.306 1.072246e-19
We explore enrichment in gene sets for HS population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC > 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in IFE granular compared to IFE spinous")
list_results[[group_name]]$enrichr_hs = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for HD population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC < 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Down-regulated in IFE granular compared to IFE spinous")
list_results[[group_name]]$enrichr_hd = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(sobj, assay = "RNA", slot = "counts"),
group1 = colnames(sobj)[sobj@active.ident %in% "IFE granular"],
group2 = colnames(sobj)[sobj@active.ident %in% "IFE spinous"])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We compute a fold change table to make a wordcloud :
fold_change = gene_corresp
colnames(fold_change) = c("gene_name", "ID")
rownames(fold_change) = fold_change$ID
fold_change$FC = ranked_gene_list[rownames(fold_change)]
fold_change$pct.1 = Seurat::GetAssayData(sobj, assay = "RNA",
slot = "counts")[, sobj@active.ident == "IFE granular"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$pct.2 = Seurat::GetAssayData(sobj, assay = "RNA",
slot = "counts")[, sobj@active.ident == "IFE spinous"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$FC_x_pct = ifelse(fold_change$FC > 0,
yes = fold_change$FC * fold_change$pct.1,
no = fold_change$FC * fold_change$pct.2)
dim(fold_change) ; head(fold_change)
## [1] 15051 6
## gene_name ID FC pct.1 pct.2
## ENSG00000238009 AL627309.1 ENSG00000238009 -1.7523930 0.000000000 0.008571429
## ENSG00000237491 AL669831.5 ENSG00000237491 -0.4392351 0.146153846 0.145714286
## ENSG00000225880 LINC00115 ENSG00000225880 -0.4404490 0.057692308 0.077142857
## ENSG00000230368 FAM41C ENSG00000230368 0.4820722 0.065384615 0.042857143
## ENSG00000230699 AL645608.3 ENSG00000230699 0.2476070 0.019230769 0.014285714
## ENSG00000187634 SAMD11 ENSG00000187634 0.2476070 0.003846154 0.002857143
## FC_x_pct
## ENSG00000238009 -0.0150205116
## ENSG00000237491 -0.0640028336
## ENSG00000225880 -0.0339774952
## ENSG00000230368 0.0315201077
## ENSG00000230699 0.0047616727
## ENSG00000187634 0.0009523345
We make the gsea plot for some gene sets :
gs_oi = c("GOBP_KERATINIZATION",
"GOCC_CORNIFIED_ENVELOPE",
"GOBP_CELLULAR_RESPONSE_TO_COPPER_ION",
"GOBP_CELLULAR_RESPONSE_TO_CADMIUM_ION",
"REACTOME_RESPONSE_TO_METAL_IONS",
"WP_TGFBETA_RECEPTOR_SIGNALING")
plot_list = make_gsea_plot(gsea_results, gs_oi, fold_change, "FC_x_pct")
patchwork::wrap_plots(plot_list, ncol = 2, widths = c(2,1.5))
We attribute score for each gene sets to all cells, to visualize the scores as violin plots.
for (one_gene_set in gs_oi) {
new_column = paste0("score_", one_gene_set)
# gene set content
gs_content = gene_sets %>%
dplyr::filter(gs_name == one_gene_set) %>%
dplyr::pull(gene_symbol) %>%
unlist()
# add score
sobj@meta.data[, new_column] = Seurat::AddModuleScore(sobj,
features = list(gs_content))$Cluster1
}
head(sobj@meta.data)
## orig.ident nCount_RNA nFeature_RNA log_nCount_RNA
## 2021_31_AAAGTGATCCTGTTGC-1 2021_31 2505 899 7.826443
## 2021_31_AACCACAAGGGCGAGA-1 2021_31 36505 5336 10.505615
## 2021_31_AAGCGAGTCGACATAC-1 2021_31 38606 5713 10.561629
## 2021_31_AATGCCAGTATCGAAA-1 2021_31 21053 4605 9.955843
## 2021_31_ACACGCGTCTAGTACG-1 2021_31 6138 1966 8.722580
## 2021_31_ACAGAAAGTTCTCCCA-1 2021_31 31947 5093 10.372240
## project_name sample_identifier sample_type
## 2021_31_AAAGTGATCCTGTTGC-1 2021_31 HS_1 HS
## 2021_31_AACCACAAGGGCGAGA-1 2021_31 HS_1 HS
## 2021_31_AAGCGAGTCGACATAC-1 2021_31 HS_1 HS
## 2021_31_AATGCCAGTATCGAAA-1 2021_31 HS_1 HS
## 2021_31_ACACGCGTCTAGTACG-1 2021_31 HS_1 HS
## 2021_31_ACAGAAAGTTCTCCCA-1 2021_31 HS_1 HS
## laboratory location Seurat.Phase cyclone.Phase
## 2021_31_AAAGTGATCCTGTTGC-1 Our axillary G2M G1
## 2021_31_AACCACAAGGGCGAGA-1 Our axillary S G1
## 2021_31_AAGCGAGTCGACATAC-1 Our axillary S G1
## 2021_31_AATGCCAGTATCGAAA-1 Our axillary G1 S
## 2021_31_ACACGCGTCTAGTACG-1 Our axillary G1 G1
## 2021_31_ACAGAAAGTTCTCCCA-1 Our axillary G1 G1
## percent.mt percent.rb cell_type
## 2021_31_AAAGTGATCCTGTTGC-1 10.774142 30.88587 IFE basal
## 2021_31_AACCACAAGGGCGAGA-1 4.364732 35.56955 IFE basal
## 2021_31_AAGCGAGTCGACATAC-1 4.020816 29.88815 sebocytes
## 2021_31_AATGCCAGTATCGAAA-1 4.317912 25.80783 IFE basal
## 2021_31_ACACGCGTCTAGTACG-1 0.228013 23.87622 IFE granular spinous
## 2021_31_ACAGAAAGTTCTCCCA-1 4.152065 29.18023 IFE granular spinous
## RNA_snn_res.0.9 seurat_clusters cluster_type
## 2021_31_AAAGTGATCCTGTTGC-1 7 7 others
## 2021_31_AACCACAAGGGCGAGA-1 6 6 IFE proliferative
## 2021_31_AAGCGAGTCGACATAC-1 0 0 IFE granular
## 2021_31_AATGCCAGTATCGAAA-1 0 0 IFE granular
## 2021_31_ACACGCGTCTAGTACG-1 4 4 IFE granular
## 2021_31_ACAGAAAGTTCTCCCA-1 2 2 IFE spinous
## score_GOBP_KERATINIZATION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.13545620
## 2021_31_AACCACAAGGGCGAGA-1 -0.05531981
## 2021_31_AAGCGAGTCGACATAC-1 -0.04239343
## 2021_31_AATGCCAGTATCGAAA-1 -0.03069908
## 2021_31_ACACGCGTCTAGTACG-1 0.03061087
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.02291899
## score_GOCC_CORNIFIED_ENVELOPE
## 2021_31_AAAGTGATCCTGTTGC-1 -0.12450299
## 2021_31_AACCACAAGGGCGAGA-1 -0.14821611
## 2021_31_AAGCGAGTCGACATAC-1 -0.19455109
## 2021_31_AATGCCAGTATCGAAA-1 0.05711020
## 2021_31_ACACGCGTCTAGTACG-1 0.49044861
## 2021_31_ACAGAAAGTTCTCCCA-1 0.02982518
## score_GOBP_CELLULAR_RESPONSE_TO_COPPER_ION
## 2021_31_AAAGTGATCCTGTTGC-1 0.15733893
## 2021_31_AACCACAAGGGCGAGA-1 -0.03086286
## 2021_31_AAGCGAGTCGACATAC-1 -0.10995941
## 2021_31_AATGCCAGTATCGAAA-1 0.01772833
## 2021_31_ACACGCGTCTAGTACG-1 -0.06047799
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.06338896
## score_GOBP_CELLULAR_RESPONSE_TO_CADMIUM_ION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.013775233
## 2021_31_AACCACAAGGGCGAGA-1 0.072102378
## 2021_31_AAGCGAGTCGACATAC-1 -0.060180186
## 2021_31_AATGCCAGTATCGAAA-1 0.111538469
## 2021_31_ACACGCGTCTAGTACG-1 -0.078435331
## 2021_31_ACAGAAAGTTCTCCCA-1 0.007335939
## score_REACTOME_RESPONSE_TO_METAL_IONS
## 2021_31_AAAGTGATCCTGTTGC-1 0.17842370
## 2021_31_AACCACAAGGGCGAGA-1 -0.11334618
## 2021_31_AAGCGAGTCGACATAC-1 -0.12670757
## 2021_31_AATGCCAGTATCGAAA-1 0.03215326
## 2021_31_ACACGCGTCTAGTACG-1 -0.20175066
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.08283927
## score_WP_TGFBETA_RECEPTOR_SIGNALING
## 2021_31_AAAGTGATCCTGTTGC-1 -0.038972501
## 2021_31_AACCACAAGGGCGAGA-1 0.038135302
## 2021_31_AAGCGAGTCGACATAC-1 -0.047914885
## 2021_31_AATGCCAGTATCGAAA-1 0.056690016
## 2021_31_ACACGCGTCTAGTACG-1 0.005909078
## 2021_31_ACAGAAAGTTCTCCCA-1 0.018563038
We visualize the scores as violin plots :
plot_list = lapply(gs_oi, FUN = function(one_gs) {
p = Seurat::VlnPlot(sobj, features = paste0("score_", one_gs),
group.by = "cluster_type",
split.by = "sample_type") +
ggplot2::scale_fill_manual(values = c("#C55F40", "#2C78E6"),
breaks = c("HS", "HD")) +
ggplot2::labs(title = one_gs) +
ggplot2::theme(axis.title.x = element_blank())
return(p)
})
patchwork::wrap_plots(plot_list, ncol = 2)
We perform a differential expression analysis for all clusters, between HS patients (HS) and healthy donors (HD).
group_name = "HS_vs_HD"
Seurat::Idents(sobj) = sobj$sample_type
table(Seurat::Idents(sobj))
##
## HS HD
## 857 137
We identify DE genes between HS and HD :
mark = Seurat::FindMarkers(sobj, ident.1 = "HS", ident.2 = "HD")
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 125 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## S100A7 1.707978e-26 1.5712922 0.797 0.314 2.570678e-22
## LGALS7B 2.619916e-20 1.4849381 0.875 0.774 3.943235e-16
## LGALS7 5.644298e-20 1.4848100 0.587 0.219 8.495234e-16
## S100A9 2.213672e-22 1.2245197 0.883 0.562 3.331798e-18
## SERPINB4 1.691713e-09 1.2032478 0.263 0.029 2.546198e-05
## KRT16 3.204700e-09 1.0306455 0.518 0.219 4.823395e-05
## MIF 5.485978e-23 0.8920172 0.667 0.321 8.256945e-19
## FABP5 8.993906e-15 0.8903571 0.918 0.759 1.353673e-10
## RPS26 3.180425e-47 0.8234142 0.982 0.993 4.786857e-43
## MTRNR2L8 4.560037e-16 0.8078924 0.840 0.774 6.863312e-12
## AKR1B10 6.070274e-28 0.7973171 0.569 0.036 9.136369e-24
## SERPINB3 2.661133e-08 0.7824037 0.219 0.015 4.005271e-04
## S100A8 1.038251e-16 0.7760133 0.887 0.613 1.562672e-12
## CRABP2 2.107146e-09 0.5757340 0.659 0.409 3.171465e-05
## MTRNR2L10 5.297609e-10 0.5584586 0.627 0.460 7.973432e-06
## HSPA8 2.017428e-14 0.4930622 0.882 0.752 3.036431e-10
## ARF5 1.640443e-21 0.4903324 0.516 0.080 2.469031e-17
## IFITM3 1.908740e-06 0.4180892 0.769 0.620 2.872845e-02
## MTRNR2L12 8.389870e-07 0.4153288 0.870 0.781 1.262759e-02
## S100A11 2.884911e-17 0.4150983 0.984 0.978 4.342079e-13
We explore enrichment in gene sets for HS population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC > 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hs = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for HD population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC < 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Down-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hd = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(sobj, assay = "RNA", slot = "counts"),
group1 = colnames(sobj)[sobj@active.ident %in% "HS"],
group2 = colnames(sobj)[sobj@active.ident %in% "HD"])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We compute a fold change table to make a wordcloud :
fold_change = gene_corresp
colnames(fold_change) = c("gene_name", "ID")
rownames(fold_change) = fold_change$ID
fold_change$FC = ranked_gene_list[rownames(fold_change)]
fold_change$pct.1 = Seurat::GetAssayData(sobj, assay = "RNA", slot = "counts")[, sobj@active.ident == "HS"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$pct.2 = Seurat::GetAssayData(sobj, assay = "RNA", slot = "counts")[, sobj@active.ident == "HD"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$FC_x_pct = ifelse(fold_change$FC > 0,
yes = fold_change$FC * fold_change$pct.1,
no = fold_change$FC * fold_change$pct.2)
dim(fold_change) ; head(fold_change)
## [1] 15051 6
## gene_name ID FC pct.1 pct.2
## ENSG00000238009 AL627309.1 ENSG00000238009 0.03273644 0.010501750 0.00000000
## ENSG00000237491 AL669831.5 ENSG00000237491 0.17957783 0.149358226 0.08029197
## ENSG00000225880 LINC00115 ENSG00000225880 0.47634309 0.071178530 0.02919708
## ENSG00000230368 FAM41C ENSG00000230368 1.14821366 0.066511085 0.01459854
## ENSG00000230699 AL645608.3 ENSG00000230699 -0.11926665 0.019836639 0.00729927
## ENSG00000187634 SAMD11 ENSG00000187634 -3.45911666 0.005834306 0.03649635
## FC_x_pct
## ENSG00000238009 0.0003437899
## ENSG00000237491 0.0268214259
## ENSG00000225880 0.0339054009
## ENSG00000230368 0.0763689364
## ENSG00000230699 -0.0008705595
## ENSG00000187634 -0.1262451334
We make the gsea plot for some gene sets :
gs_oi = c("HALLMARK_INTERFERON_ALPHA_RESPONSE",
"HALLMARK_INTERFERON_GAMMA_RESPONSE",
"REACTOME_KERATINIZATION",
"GOBP_KERATINIZATION",
"GOBP_MACROPHAGE_CYTOKINE_PRODUCTION",
"GOCC_MHC_PROTEIN_COMPLEX",
"KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION")
plot_list = make_gsea_plot(gsea_results, gs_oi, fold_change, "FC_x_pct")
patchwork::wrap_plots(plot_list, ncol = 2, widths = c(2,1.5))
We attribute score for each gene sets to all cells, to visualize the scores as violin plots.
for (one_gene_set in gs_oi) {
new_column = paste0("score_", one_gene_set)
# gene set content
gs_content = gene_sets %>%
dplyr::filter(gs_name == one_gene_set) %>%
dplyr::pull(gene_symbol) %>%
unlist()
# add score
sobj@meta.data[, new_column] = Seurat::AddModuleScore(sobj,
features = list(gs_content))$Cluster1
}
head(sobj@meta.data)
## orig.ident nCount_RNA nFeature_RNA log_nCount_RNA
## 2021_31_AAAGTGATCCTGTTGC-1 2021_31 2505 899 7.826443
## 2021_31_AACCACAAGGGCGAGA-1 2021_31 36505 5336 10.505615
## 2021_31_AAGCGAGTCGACATAC-1 2021_31 38606 5713 10.561629
## 2021_31_AATGCCAGTATCGAAA-1 2021_31 21053 4605 9.955843
## 2021_31_ACACGCGTCTAGTACG-1 2021_31 6138 1966 8.722580
## 2021_31_ACAGAAAGTTCTCCCA-1 2021_31 31947 5093 10.372240
## project_name sample_identifier sample_type
## 2021_31_AAAGTGATCCTGTTGC-1 2021_31 HS_1 HS
## 2021_31_AACCACAAGGGCGAGA-1 2021_31 HS_1 HS
## 2021_31_AAGCGAGTCGACATAC-1 2021_31 HS_1 HS
## 2021_31_AATGCCAGTATCGAAA-1 2021_31 HS_1 HS
## 2021_31_ACACGCGTCTAGTACG-1 2021_31 HS_1 HS
## 2021_31_ACAGAAAGTTCTCCCA-1 2021_31 HS_1 HS
## laboratory location Seurat.Phase cyclone.Phase
## 2021_31_AAAGTGATCCTGTTGC-1 Our axillary G2M G1
## 2021_31_AACCACAAGGGCGAGA-1 Our axillary S G1
## 2021_31_AAGCGAGTCGACATAC-1 Our axillary S G1
## 2021_31_AATGCCAGTATCGAAA-1 Our axillary G1 S
## 2021_31_ACACGCGTCTAGTACG-1 Our axillary G1 G1
## 2021_31_ACAGAAAGTTCTCCCA-1 Our axillary G1 G1
## percent.mt percent.rb cell_type
## 2021_31_AAAGTGATCCTGTTGC-1 10.774142 30.88587 IFE basal
## 2021_31_AACCACAAGGGCGAGA-1 4.364732 35.56955 IFE basal
## 2021_31_AAGCGAGTCGACATAC-1 4.020816 29.88815 sebocytes
## 2021_31_AATGCCAGTATCGAAA-1 4.317912 25.80783 IFE basal
## 2021_31_ACACGCGTCTAGTACG-1 0.228013 23.87622 IFE granular spinous
## 2021_31_ACAGAAAGTTCTCCCA-1 4.152065 29.18023 IFE granular spinous
## RNA_snn_res.0.9 seurat_clusters cluster_type
## 2021_31_AAAGTGATCCTGTTGC-1 7 7 others
## 2021_31_AACCACAAGGGCGAGA-1 6 6 IFE proliferative
## 2021_31_AAGCGAGTCGACATAC-1 0 0 IFE granular
## 2021_31_AATGCCAGTATCGAAA-1 0 0 IFE granular
## 2021_31_ACACGCGTCTAGTACG-1 4 4 IFE granular
## 2021_31_ACAGAAAGTTCTCCCA-1 2 2 IFE spinous
## score_GOBP_KERATINIZATION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.13545620
## 2021_31_AACCACAAGGGCGAGA-1 -0.05531981
## 2021_31_AAGCGAGTCGACATAC-1 -0.04239343
## 2021_31_AATGCCAGTATCGAAA-1 -0.03069908
## 2021_31_ACACGCGTCTAGTACG-1 0.03061087
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.02291899
## score_GOCC_CORNIFIED_ENVELOPE
## 2021_31_AAAGTGATCCTGTTGC-1 -0.12450299
## 2021_31_AACCACAAGGGCGAGA-1 -0.14821611
## 2021_31_AAGCGAGTCGACATAC-1 -0.19455109
## 2021_31_AATGCCAGTATCGAAA-1 0.05711020
## 2021_31_ACACGCGTCTAGTACG-1 0.49044861
## 2021_31_ACAGAAAGTTCTCCCA-1 0.02982518
## score_GOBP_CELLULAR_RESPONSE_TO_COPPER_ION
## 2021_31_AAAGTGATCCTGTTGC-1 0.15733893
## 2021_31_AACCACAAGGGCGAGA-1 -0.03086286
## 2021_31_AAGCGAGTCGACATAC-1 -0.10995941
## 2021_31_AATGCCAGTATCGAAA-1 0.01772833
## 2021_31_ACACGCGTCTAGTACG-1 -0.06047799
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.06338896
## score_GOBP_CELLULAR_RESPONSE_TO_CADMIUM_ION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.013775233
## 2021_31_AACCACAAGGGCGAGA-1 0.072102378
## 2021_31_AAGCGAGTCGACATAC-1 -0.060180186
## 2021_31_AATGCCAGTATCGAAA-1 0.111538469
## 2021_31_ACACGCGTCTAGTACG-1 -0.078435331
## 2021_31_ACAGAAAGTTCTCCCA-1 0.007335939
## score_REACTOME_RESPONSE_TO_METAL_IONS
## 2021_31_AAAGTGATCCTGTTGC-1 0.17842370
## 2021_31_AACCACAAGGGCGAGA-1 -0.11334618
## 2021_31_AAGCGAGTCGACATAC-1 -0.12670757
## 2021_31_AATGCCAGTATCGAAA-1 0.03215326
## 2021_31_ACACGCGTCTAGTACG-1 -0.20175066
## 2021_31_ACAGAAAGTTCTCCCA-1 -0.08283927
## score_WP_TGFBETA_RECEPTOR_SIGNALING
## 2021_31_AAAGTGATCCTGTTGC-1 -0.038972501
## 2021_31_AACCACAAGGGCGAGA-1 0.038135302
## 2021_31_AAGCGAGTCGACATAC-1 -0.047914885
## 2021_31_AATGCCAGTATCGAAA-1 0.056690016
## 2021_31_ACACGCGTCTAGTACG-1 0.005909078
## 2021_31_ACAGAAAGTTCTCCCA-1 0.018563038
## score_HALLMARK_INTERFERON_ALPHA_RESPONSE
## 2021_31_AAAGTGATCCTGTTGC-1 -0.01163969
## 2021_31_AACCACAAGGGCGAGA-1 0.04052888
## 2021_31_AAGCGAGTCGACATAC-1 0.02563452
## 2021_31_AATGCCAGTATCGAAA-1 0.07830894
## 2021_31_ACACGCGTCTAGTACG-1 -0.07371936
## 2021_31_ACAGAAAGTTCTCCCA-1 0.05225323
## score_HALLMARK_INTERFERON_GAMMA_RESPONSE
## 2021_31_AAAGTGATCCTGTTGC-1 -0.03163757
## 2021_31_AACCACAAGGGCGAGA-1 0.01897797
## 2021_31_AAGCGAGTCGACATAC-1 0.01765398
## 2021_31_AATGCCAGTATCGAAA-1 0.04969244
## 2021_31_ACACGCGTCTAGTACG-1 -0.04287433
## 2021_31_ACAGAAAGTTCTCCCA-1 0.02979381
## score_REACTOME_KERATINIZATION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.004632692
## 2021_31_AACCACAAGGGCGAGA-1 0.057549695
## 2021_31_AAGCGAGTCGACATAC-1 0.066181820
## 2021_31_AATGCCAGTATCGAAA-1 0.157580089
## 2021_31_ACACGCGTCTAGTACG-1 0.407402733
## 2021_31_ACAGAAAGTTCTCCCA-1 0.198902298
## score_GOBP_MACROPHAGE_CYTOKINE_PRODUCTION
## 2021_31_AAAGTGATCCTGTTGC-1 -0.06397283
## 2021_31_AACCACAAGGGCGAGA-1 0.09441426
## 2021_31_AAGCGAGTCGACATAC-1 0.05652529
## 2021_31_AATGCCAGTATCGAAA-1 0.05181899
## 2021_31_ACACGCGTCTAGTACG-1 0.03751105
## 2021_31_ACAGAAAGTTCTCCCA-1 0.04375760
## score_GOCC_MHC_PROTEIN_COMPLEX
## 2021_31_AAAGTGATCCTGTTGC-1 -0.01578330
## 2021_31_AACCACAAGGGCGAGA-1 0.36108395
## 2021_31_AAGCGAGTCGACATAC-1 0.16194305
## 2021_31_AATGCCAGTATCGAAA-1 0.39295972
## 2021_31_ACACGCGTCTAGTACG-1 0.03782853
## 2021_31_ACAGAAAGTTCTCCCA-1 0.34138984
## score_KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION
## 2021_31_AAAGTGATCCTGTTGC-1 0.06097775
## 2021_31_AACCACAAGGGCGAGA-1 0.20508977
## 2021_31_AAGCGAGTCGACATAC-1 0.24957622
## 2021_31_AATGCCAGTATCGAAA-1 0.26696440
## 2021_31_ACACGCGTCTAGTACG-1 0.09250680
## 2021_31_ACAGAAAGTTCTCCCA-1 0.17942810
We visualize the scores as violin plots :
plot_list = lapply(gs_oi, FUN = function(one_gs) {
p = Seurat::VlnPlot(sobj, features = paste0("score_", one_gene_set),
group.by = "cluster_type",
split.by = "sample_type") +
ggplot2::scale_fill_manual(values = c("#C55F40", "#2C78E6"),
breaks = c("HS", "HD")) +
ggplot2::labs(title = one_gs) +
ggplot2::theme(axis.title.x = element_blank())
return(p)
})
patchwork::wrap_plots(plot_list, ncol = 2)
We represent differentially expressed genes. First, we extract some DE genes :
features_oi = lapply(list_results, `[[`, "mark") %>%
do.call(rbind.data.frame, .) %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::filter(pct.1 - pct.2 > 0.2 | abs(avg_logFC) > 1) %>%
rownames() %>%
stringr::str_split(., pattern = "\\.") %>%
lapply(., `[[`, 2) %>%
unlist() %>%
unique()
length(features_oi)
## [1] 95
We prepare the scaled expression matrix :
mat_expression = Seurat::GetAssayData(sobj, assay = "RNA", slot = "data")[features_oi, ]
mat_expression = Matrix::t(mat_expression)
mat_expression = dynutils::scale_quantile(mat_expression) # between 0 and 1
mat_expression = Matrix::t(mat_expression)
mat_expression = as.matrix(mat_expression) # not sparse
dim(mat_expression)
## [1] 95 994
We prepare the heatmap annotation :
ha_top = ComplexHeatmap::HeatmapAnnotation(
sample_type = sobj$sample_type,
cluster_type = sobj$cluster_type,
cluster = sobj$seurat_clusters,
col = list(sample_type = setNames(nm = c("HS", "HD"),
c("#C55F40", "#2C78E6")),
cluster_type = cluster_color,
cluster = setNames(nm = levels(sobj$seurat_clusters),
aquarius::gg_color_hue(length(levels(sobj$seurat_clusters))))))
And the heatmap :
ht = ComplexHeatmap::Heatmap(mat_expression,
col = aquarius::color_cnv,
# Annotation
top_annotation = ha_top,
# Grouping
column_order = sobj@meta.data %>%
dplyr::arrange(cluster_type, sample_type, seurat_clusters) %>%
rownames(),
column_title = NULL,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_column_names = FALSE,
# Visual aspect
show_heatmap_legend = TRUE,
border = TRUE)
ComplexHeatmap::draw(ht,
merge_legend = TRUE,
heatmap_legend_side = "bottom",
annotation_legend_side = "bottom")
We save the list of results :
saveRDS(list_results, file = paste0(out_dir, "/", save_name, "_list_results.rds"))
We save the Seurat object with the new annotation :
saveRDS(sobj, file = paste0(out_dir, "/", save_name, "_sobj_annotated.rds"))
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## locale:
## [1] C
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ComplexHeatmap_2.14.0 ggplot2_3.3.5 patchwork_1.1.2
## [4] dplyr_1.0.7
##
## loaded via a namespace (and not attached):
## [1] softImpute_1.4 graphlayouts_0.7.0
## [3] pbapply_1.4-2 lattice_0.20-41
## [5] haven_2.3.1 vctrs_0.3.8
## [7] usethis_2.0.1 dynwrap_1.2.1
## [9] blob_1.2.1 survival_3.2-13
## [11] prodlim_2019.11.13 dynutils_1.0.5
## [13] later_1.3.0 DBI_1.1.1
## [15] R.utils_2.11.0 SingleCellExperiment_1.8.0
## [17] rappdirs_0.3.3 uwot_0.1.8
## [19] dqrng_0.2.1 jpeg_0.1-8.1
## [21] zlibbioc_1.32.0 pspline_1.0-18
## [23] pcaMethods_1.78.0 mvtnorm_1.1-1
## [25] htmlwidgets_1.5.4 GlobalOptions_0.1.2
## [27] future_1.22.1 UpSetR_1.4.0
## [29] laeken_0.5.2 leiden_0.3.3
## [31] clustree_0.4.3 parallel_3.6.3
## [33] scater_1.14.6 irlba_2.3.3
## [35] markdown_1.1 DEoptimR_1.0-9
## [37] tidygraph_1.1.2 Rcpp_1.0.9
## [39] readr_2.0.2 KernSmooth_2.23-17
## [41] carrier_0.1.0 promises_1.1.0
## [43] gdata_2.18.0 DelayedArray_0.12.3
## [45] limma_3.42.2 graph_1.64.0
## [47] RcppParallel_5.1.4 Hmisc_4.4-0
## [49] fs_1.5.2 RSpectra_0.16-0
## [51] fastmatch_1.1-0 ranger_0.12.1
## [53] digest_0.6.25 png_0.1-7
## [55] sctransform_0.2.1 cowplot_1.0.0
## [57] DOSE_3.12.0 here_1.0.1
## [59] TInGa_0.0.0.9000 ggraph_2.0.3
## [61] pkgconfig_2.0.3 GO.db_3.10.0
## [63] DelayedMatrixStats_1.8.0 gower_0.2.1
## [65] ggbeeswarm_0.6.0 iterators_1.0.12
## [67] DropletUtils_1.6.1 reticulate_1.26
## [69] clusterProfiler_3.14.3 SummarizedExperiment_1.16.1
## [71] circlize_0.4.15 beeswarm_0.4.0
## [73] GetoptLong_1.0.5 xfun_0.35
## [75] bslib_0.3.1 zoo_1.8-10
## [77] tidyselect_1.1.0 reshape2_1.4.4
## [79] purrr_0.3.4 ica_1.0-2
## [81] pcaPP_1.9-73 viridisLite_0.3.0
## [83] rtracklayer_1.46.0 rlang_1.0.2
## [85] hexbin_1.28.1 jquerylib_0.1.4
## [87] dyneval_0.9.9 glue_1.4.2
## [89] RColorBrewer_1.1-2 matrixStats_0.56.0
## [91] stringr_1.4.0 lava_1.6.7
## [93] europepmc_0.3 DESeq2_1.26.0
## [95] recipes_0.1.17 labeling_0.3
## [97] httpuv_1.5.2 class_7.3-17
## [99] BiocNeighbors_1.4.2 DO.db_2.9
## [101] annotate_1.64.0 jsonlite_1.7.2
## [103] XVector_0.26.0 bit_4.0.4
## [105] mime_0.9 aquarius_0.1.5
## [107] Rsamtools_2.2.3 gridExtra_2.3
## [109] gplots_3.0.3 stringi_1.4.6
## [111] processx_3.5.2 gsl_2.1-6
## [113] bitops_1.0-6 cli_3.0.1
## [115] batchelor_1.2.4 RSQLite_2.2.0
## [117] randomForest_4.6-14 tidyr_1.1.4
## [119] data.table_1.14.2 rstudioapi_0.13
## [121] org.Mm.eg.db_3.10.0 GenomicAlignments_1.22.1
## [123] nlme_3.1-147 qvalue_2.18.0
## [125] scran_1.14.6 locfit_1.5-9.4
## [127] scDblFinder_1.1.8 listenv_0.8.0
## [129] ggthemes_4.2.4 gridGraphics_0.5-0
## [131] R.oo_1.24.0 dbplyr_1.4.4
## [133] BiocGenerics_0.32.0 TTR_0.24.2
## [135] readxl_1.3.1 lifecycle_1.0.1
## [137] timeDate_3043.102 ggpattern_0.3.1
## [139] munsell_0.5.0 cellranger_1.1.0
## [141] R.methodsS3_1.8.1 proxyC_0.1.5
## [143] visNetwork_2.0.9 caTools_1.18.0
## [145] codetools_0.2-16 ggwordcloud_0.5.0
## [147] Biobase_2.46.0 GenomeInfoDb_1.22.1
## [149] vipor_0.4.5 lmtest_0.9-38
## [151] msigdbr_7.5.1 htmlTable_1.13.3
## [153] triebeard_0.3.0 lsei_1.2-0
## [155] xtable_1.8-4 ROCR_1.0-7
## [157] BiocManager_1.30.10 scatterplot3d_0.3-41
## [159] abind_1.4-5 farver_2.0.3
## [161] parallelly_1.28.1 RANN_2.6.1
## [163] askpass_1.1 GenomicRanges_1.38.0
## [165] RcppAnnoy_0.0.16 tibble_3.1.5
## [167] ggdendro_0.1-20 cluster_2.1.0
## [169] future.apply_1.5.0 Seurat_3.1.5
## [171] dendextend_1.15.1 Matrix_1.3-2
## [173] ellipsis_0.3.2 prettyunits_1.1.1
## [175] lubridate_1.7.9 ggridges_0.5.2
## [177] igraph_1.2.5 RcppEigen_0.3.3.7.0
## [179] fgsea_1.12.0 remotes_2.4.2
## [181] scBFA_1.0.0 destiny_3.0.1
## [183] VIM_6.1.1 testthat_3.1.0
## [185] htmltools_0.5.2 BiocFileCache_1.10.2
## [187] yaml_2.2.1 utf8_1.1.4
## [189] plotly_4.9.2.1 XML_3.99-0.3
## [191] ModelMetrics_1.2.2.2 e1071_1.7-3
## [193] foreign_0.8-76 withr_2.5.0
## [195] fitdistrplus_1.0-14 BiocParallel_1.20.1
## [197] xgboost_1.4.1.1 bit64_4.0.5
## [199] foreach_1.5.0 robustbase_0.93-9
## [201] Biostrings_2.54.0 GOSemSim_2.13.1
## [203] rsvd_1.0.3 memoise_2.0.0
## [205] evaluate_0.18 forcats_0.5.0
## [207] rio_0.5.16 geneplotter_1.64.0
## [209] tzdb_0.1.2 caret_6.0-86
## [211] ps_1.6.0 DiagrammeR_1.0.6.1
## [213] curl_4.3 fdrtool_1.2.15
## [215] fansi_0.4.1 highr_0.8
## [217] urltools_1.7.3 xts_0.12.1
## [219] GSEABase_1.48.0 acepack_1.4.1
## [221] edgeR_3.28.1 checkmate_2.0.0
## [223] scds_1.2.0 cachem_1.0.6
## [225] npsurv_0.4-0 babelgene_22.3
## [227] rjson_0.2.20 openxlsx_4.1.5
## [229] ggrepel_0.9.1 clue_0.3-60
## [231] rprojroot_2.0.2 stabledist_0.7-1
## [233] tools_3.6.3 sass_0.4.0
## [235] nichenetr_1.1.1 magrittr_2.0.1
## [237] RCurl_1.98-1.2 proxy_0.4-24
## [239] car_3.0-11 ape_5.3
## [241] ggplotify_0.0.5 xml2_1.3.2
## [243] httr_1.4.2 assertthat_0.2.1
## [245] rmarkdown_2.18 boot_1.3-25
## [247] globals_0.14.0 R6_2.4.1
## [249] Rhdf5lib_1.8.0 nnet_7.3-14
## [251] RcppHNSW_0.2.0 progress_1.2.2
## [253] genefilter_1.68.0 statmod_1.4.34
## [255] gtools_3.8.2 shape_1.4.6
## [257] HDF5Array_1.14.4 BiocSingular_1.2.2
## [259] rhdf5_2.30.1 splines_3.6.3
## [261] AUCell_1.8.0 carData_3.0-4
## [263] colorspace_1.4-1 generics_0.1.0
## [265] stats4_3.6.3 base64enc_0.1-3
## [267] dynfeature_1.0.0 smoother_1.1
## [269] gridtext_0.1.1 pillar_1.6.3
## [271] tweenr_1.0.1 sp_1.4-1
## [273] ggplot.multistats_1.0.0 rvcheck_0.1.8
## [275] GenomeInfoDbData_1.2.2 plyr_1.8.6
## [277] gtable_0.3.0 zip_2.2.0
## [279] knitr_1.41 latticeExtra_0.6-29
## [281] biomaRt_2.42.1 IRanges_2.20.2
## [283] fastmap_1.1.0 ADGofTest_0.3
## [285] copula_1.0-0 doParallel_1.0.15
## [287] AnnotationDbi_1.48.0 vcd_1.4-8
## [289] babelwhale_1.0.1 openssl_1.4.1
## [291] scales_1.1.1 backports_1.2.1
## [293] S4Vectors_0.24.4 ipred_0.9-12
## [295] enrichplot_1.6.1 hms_1.1.1
## [297] ggforce_0.3.1 Rtsne_0.15
## [299] shiny_1.7.1 numDeriv_2016.8-1.1
## [301] polyclip_1.10-0 lazyeval_0.2.2
## [303] Formula_1.2-3 tsne_0.1-3
## [305] crayon_1.3.4 MASS_7.3-54
## [307] pROC_1.16.2 viridis_0.5.1
## [309] dynparam_1.0.0 rpart_4.1-15
## [311] zinbwave_1.8.0 compiler_3.6.3
## [313] ggtext_0.1.0